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Abstract — The linear noise approximation is commonly used 
to obtain intrinsic noise statistics such as Fano factors and coef- 
ficients of variation for biochemical networks. These estimates 
are accurate for networks with large numbers of molecules. 
However it is well known that many biochemical networks 
are characterized by at least one species with a small number 
of molecules. We here describe modifications to the software 
intrinsic Noise Analyzer (iNA) which enable it to accurately 
compute noise statistics over wide ranges of molecule numbers. 
This is achieved by calculating the next order corrections to 
the linear noise approximation's estimates of variance and 
covariance of concentration fluctuations. The efficiency of 
the methods is significantly improved by automated just-in- 
time compilation using the LLVM framework leading to a 
fluctuation analysis which typically outperforms that obtained 
by means of exact stochastic simulations. iNA is hence partic- 
ularly wefl suited for the needs of the computational biology 
community. 

KeywordsSiQch^siic modeling; Linear Noise Approximation; 
genetic regulatory circuits 

I. Introduction 

Experimental studies have shown that the protein abun- 
dance varies from few tens to several thousands per protein 
species per cell |1|. It is also known that the standard 
deviation of the concentration fluctuations due to the random 
timing of molecular events (intrinsic noise) roughly scales as 
the square root of the mean number of molecules |2|. Hence 
it is expected that intrinsic noise plays an important role in 
the dynamics of those biochemical networks characterized 
by at least one species with low molecule numbers. 

The stochastic simulation algorithm (SSA) is the conven- 
tional means of probing stochasticity in biochemical reaction 
systems |[3| . This method simulates every reaction event and 
hence is typically slow for large reaction networks; this 
is particularly true if one is interested in intrinsic noise 
statistics which require considerable ensemble averaging of 
the trajectories produced by the SSA. A different route of 
inferring the required statistics involves finding an approxi- 
mate solution of the chemical master equation (CME), a set 
of differential equations for the probabilities of the states of 
the system, which is mathematically equivalent to the SSA. 



We recently developed intrinsic Noise Analyzer (iNA) ||4|, 
the first software package enabling a fluctuation analysis of 
biochemical networks via the Linear Noise Approximation 
(LNA) and Effective Mesoscopic Rate Equation (EMRE) 
approximations of the CME. The former gives the variance 
and covariance of concentration fluctuations in the limit of 
large molecule numbers while the latter gives the mean 
concentrations for intermediate to large molecule numbers 
and is more accurate than the conventional Rate Equations 
(REs). 

In this paper we develop and efficiently implement in 
iNA, the Inverse Omega Square (lOS) method which gives 
accurate estimates for the mean concentrations and vari- 
ance / covariance of noise about them for systems whose 
molecule numbers vary over wide ranges (few to thousands 
of molecules). The new method is tested on two gene expres- 
sion models involving bimolecular reactions. Remarkably, 
the results of a few seconds long lOS calculation is found 
to agree very well with those from hour long ensemble 
averaging of SSA trajectories. 

II. Results 

In this section we describe the results of the novel lOS 
method implemented inside iNA, compare with the results 
of the RE, LNA and EMRE approximations of the CME and 
with exact stochastic simulations using the SSA and finally 
discuss the implementation and its performance. The three 
methods (LNA, EMRE, lOS) are derived from the system- 
size expansion (SSE) of the CME |2| which is applicable 
to monostable chemical systems. Technical details of the 
various approximation methods are provided in the section 
Methods. 

A. Applications 

We consider a two-stage model of gene expression with 
an enzymatic degradation mechanism: 

Gene ^ Gene + mRNA, mRNA ^ 0, 
mRNA ^ mRNA + Protein, 
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Table I: Kinetic parameters used for the gene expression 
model, scheme [T] as discussed in the main text. The volume 
is fixed to ft = 10~^^/ with an enzyme concentration 
of 0.1 /iM corresponding to 60 enzyme molecules. The 
Michaelis-Menten constant is 0.01/iM in all cases. 



ki 

Protein + Enzyme Complex — > Enzyme + 0. (1) 

k-i 

The scheme involves the transcription of mRNA, its trans- 
lation to protein and subsequent degradation of both mRNA 
and protein. Note that while mRNA is degraded via an un- 
specific linear reaction, the degradation of protein occurs via 
an enzyme catalyzed reaction. The latter may model proteol- 
ysis, the consumption of protein by a metabolic pathway or 
other post-translational modifications. A simplified version 
of this model is one in which the protein degradation is 
replaced by the linear reaction: Protein 0. Over the past 
decade the latter model has been the subject of numerous 
studies, principally because it can be solved exactly since the 
scheme is composed of purely first-order reactions |[5|-||7|. 
However, the former model as given by scheme ([T]), cannot 
be solved exactly because of the bimolecular association 
reaction between enzyme and protein. Hence in what follows 
we demonstrate the power of approximation methods to infer 
useful information regarding the mechanism's intrinsic noise 
properties. 

We consider the model with three different parameter sets 
(see Tab|l]) and fix the compartment volume to one femtoliter 
(one micron cubed). For all three cases, the REs predict 
the same steady state mRNA and protein concentrations: 
[mRNA] = 0.12/iM and [Protein] = 0.09/iM. These 
correspond to 72 and 54 molecule numbers, respectively. 

We have used iNA to compute the mean concentrations 
using the REs and the variance of fluctuations using the 
LNA for parameter set (i) in which transcription is fast. 
Comparing these with SSA estimates (see Fig. [T]) we see 
that the REs and LNA provide reasonably accurate results 
for this parameter set. This analysis was within the scope 
of the previous version of iNA |4|. However, the scenario 
considered is not particularly realistic. This is since the 
ratio of protein and mRNA lifetimes in this example is 
approximately 100 (as estimated from the time taken for 
the concentrations to reach 90% of their steady state values) 
while an evaluation of 1 ,962 genes in budding yeast showed 
that the ratios have median and mode close to 3 |7 |. 

We now consider parameter set (ii), the case of moderate 
transcription. In Fig. [2ja) and (b) we compare the RE and 
LNA predictions of mean concentrations and variance of 



fluctuations with that obtained from the SSA. Notice that 
in this case, the two are in severe disagreement. The SSA 
predicts that the mean concentration of protein is larger than 
that of the mRNA while the REs predict the opposite. It 
is also the case that the variance estimate of the LNA is 
considerably smaller than that of the SSA. In Fig. [2jc) we 
show the mean concentrations computed using the EMRE 
and the variance computed using the lOS method. Note that 
the latter are in good agreement with the SSA in Fig. |2jb). 
Note further that while the EMRE was already implemented 
in the previous version of iNA, the lOS was not. Hence 
the present version is the first to provide estimates of the 
mean concentrations and of the size of the noise which 
are more accurate than both the REs and the LNA. The 
transient times for this case are given by 37 minutes for 
protein and 12 minutes for mRNA concentrations with a 
ratio of approximately 3 in agreement with the median and 
mode of experimentally measured ratios. Hence this example 
provides clear evidence of the need to go beyond the RE 
and LNA level of approximation for physiologically relevant 
parameters of the gene expression model. 

Finally we consider parameter set (iii), the case of moder- 
ate transcription and fast translation (kg is an order of mag- 
nitude larger in (iii) compared to (ii)). Previous studies have 
shown that increased translation efficiency leads to increased 
noise in the protein abundance | 5 1 due to the proteins being 
produced in bursts |7 |. Indeed a single realization of the 
SSA shows proteins expressed in sharply peaked bursts (Fig. 
[3ja)). We used iNA to compute the mean concentrations 
of mRNA and protein according to RE, EMRE and lOS 
approximation methods (Fig. [3jb)). These are contrasted 
with the same mean concentrations computed from SSA 
simulations (Fig. [Sjc)). Note that lOS provides the most 
accurate result, followed by the EMRE and the REs. The 
latter performs poorly because it does not take into account 
the effect of large fluctuations due to bursty behavior on 
the mean concentrations. Transient times of 48 minutes and 
12 minutes for protein and mRNA concentrations have been 
extracted from the time course data as for previous cases; 
these are similar to those observed in the expression of the 
E. coli proteome |8|, once again showing the necessity of 
approximation methods beyond the RE and LNA to study 
physiologically relevant cases. 

B. Implementation 

iNA's framework consist of three layers of abstraction: the 
SBML parser which sets up a mathematical representation 
of the reaction network, a module which performs the SSE 
analytically using the computer algebra system Ginac J91 
and a just-in-time (JIT) compiler based on the LLVM flOl 
framework which compiles the mathematical model into 
platform specific machine code at runtime yielding high 
performance of SSE and SSA analyses implemented in iNA. 

As elaborated in the Methods section, in addition to the 
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Figure 1: Gene expression model with fast transcription rate. In panels (a) and (b) we compare the RE predictions of mean 
concentrations with those obtained from ensemble averaging 3, 000 SSA trajectories. The shaded areas denote the region of 
one standard deviation around the average concentrations which in (a) has been computed using the LNA and in (b) using 
the SSA. The results in (a) and (b) are in approximate agreement. 
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Figure 2: Gene expression model with moderate transcription rate. In panels (a) and (b) we compare the RE and LNA 
predictions of mean concentrations and variance of fluctuations with those obtained from ensemble averaging 30,000 
stochastic realizations computed using the SSA. Note that the RE and LNA predictions are very different than the actual 
values. In panel (c) we show the mean concentration prediction according to the EMRE and the variance prediction according 
to lOS. These are in good agreement with those obtained from the SSA. 



REs, van Kampen's SSE upon which the LNA, EMRE and 
lOS methods are based, requires the numerical solution of 
the following high dimensional system of linear equations 

^ = BxssE + A, (2) 
at 

whose coefficients are automatically computed by iNA using 
the system size expansion from an SBML file. Note that 
the coefficient matrix = B_{[X]) and the vector A = 
A([X]) depend parametrically on the solution of the REs. 
The vector xsse is defined as (1^"^xlna,emre, ^"^^los) 
where xlna,emre defines the variance/covariance of the LNA 
together with the corrections to mean concentrations of the 
REs according to the EMRE and xios defines the corrections 
to the variance/covariance of the LNA and the corrections to 
the mean concentrations of the EMRE according to the lOS. 
A particular simple result can be inferred by setting the time 
derivative of the REs and Eq. ^ to zero yielding B x^^^ = 



—A whose solution is iNA's steady state analysis. First 
the typically nonlinear REs are solved using the Newton- 
Raphson method with line search |11| and consequently 
the solution to the latter linear equation is found using 
standard linear algebra. In the present version of iNA, time- 
course analysis is efficiently performed by means of the all- 
round ODE integrator LSODA which automatically switches 
between explicit and implicit methods |12|. 

The number of simultaneous equations to be solved for 
the LNA, EMRE analysis is approximately proportional 
to N"^ where N is the number of independent species 
after conservation analysis pSj while for lOS, the number 
of equations is approximately proportional to A^^. The 
quadratic and cubic dependencies represent a challenge for 
software development. iNA's previous version addressed this 
need for performance by providing automated OpenMP 
parallelism and a bytecode interpreter (BCI) for efficient 
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Figure 3: Gene expression model with efficient translation. In panel (a) we show a single SSA trajectory which illustrates 
the large protein bursts in the protein concentrations. In panel (b) we show the concentrations predicted by REs, EMRE and 
lOS methods and in (c) we show the mean concentrations obtained from ensemble averaging 30, 000 SSA trajectories. Only 
the lOS method is in good quantitative agreement with the SSA predictions. These predictions are over 5 times larger than 
those of REs, the discrepancies stemming from the fact that the latter do not take into account contributions to the mean 
due to large fluctuations. 



expression evaluation |4|. The latter concept has proven 
its performance for both SSE and SSA methods while 
maintaining compatibility over many platforms. With the 
present version we introduce a cross platform strategy for 
JIT compilation of the SSE or SSA methods. iNA uses 
the modern compiler framework LLVM providing fast JIT 
compilation that allows to emmit and execute platform 
specific machine code at runtime |10|. This allows us to 
automatically compile the system size expansion ODEs as 
native machine code executable directly on the CPU without 
resorting to interpreters. Therefore iNA's JIT feature enjoys 
the speed of statically compiled code while maintaining 
the flexibility common to interpreters. Moreover, the LLVM 
compiler framework offers the use of optimizations on 
platform independent and platform specific instruction levels 
which become advantageous for computationally expensive 
calculations. 

C. Performance 

In order to benchmark the present implementation we 
consider an ensemble of independent single cell oscillators 
with weak negative feedback involving 10 species and 16 
reactions. The model describes the transcription of mRNA, 
M, by a a clock gene, G and its translation to a cytosolic 
clock protein, Pq'. 

G%G^M, M^0, M^M + Pc. (3) 

A negative feedback loop arises from nuclear import of 
cytosolic protein and its binding to the gene which represses 
transcription: 

Pc ^ Pn, Pc ^ S, 

fcout 



Pn + G^ GPn, Pn + GPm ^ GPm2, (4) 

/c_i k-2 

where P/v denotes the the nuclear protein. Note that protein 
Pc is consumed also by a different pathway involving S 
which is not explicitly modeled. As in the preceding exam- 
ple, degradation occurs via enzyme-catalyzed mechanisms, 
in this case via cytosolic and nuclear proteases E and F: 

Pa + E^EPc^E, 

fc-3 

P^ + F^FPj,^ F. (5) 

The model has been analyzed in detail in Ref. ||4j using the 
EMREs computed by iNA. Therein it has been shown that 
intrinsic noise can amplify transient oscillations which are 
visible even at the cell population level. In Fig.|4]we compare 
the outcome of 50, 000 independent realizations of the SSA 
with iNA's prediction of these synchronous noise-induced 
oscillations using the EMRE for the mean concentrations 
together with the lOS estimate of variance (the new feature 
in the present version of iNA). These estimates agree well 
with those obtained from simulations. 

After conservation analysis the model comprises 7 species 
yielding a total number of 161 simultaneous equations for 
the SSE and hence is well suited for direct benchmarking 
purposes. This is particularly challenging for ODE integra- 
tors since the underlying stochastic dynamics causes the 
full system of 161 coupled equations to exhibit damped 
oscillations. The results of the benchmarks are summarized 
in Tab. [ll| highlighting the performance of the present version 
of iNA over the previous version. The improvements of 
iNA's SSE using the LSODA over the previous Rosenbrock 
method reduces the execution time by a factor 2 — 3 while 
factors of 5 — 6 are attained using the JIT compiler which 



combined reduce the execution time from 255 to less than 
2s. This is compared to the execution time of the SSA 
which is computationally extremely expensive because of 
the considerable number of trajectories which need to be 
averaged to obtain accurate statistics. 



method 


SSE, LSODA 


SSE, Rosenbrock 


SSA, single/ens. 


BCI 


12.88s (13.13s) 


25.40s (25.66s) 


0.15s/2.0h 


BCI (opt.) 


8.21s (8.51s) 


26.59s (26.89s) 


0.13s/1.8h 


JIT 


1.46s (1.86s) 


10.77s (15.51s) 


0.10s/1.4h 


JIT (opt.) 


1.12s (6.34s) 


10.91s (17.04s) 


0.10s/1.4h 



Table II: Execution times of iNA for SSE and SSA for 
the gene oscillatory network reproducing Fig. |4] (a) and 
(b), respectively. The table compares the SSE method using 
the ODE integrator LSODA and Rosenbrock in combination 
with a bytecode interpreter (BCI) or iNA's JIT feature. The 
value in the brackets denotes the overall time including byte- 
code assembly or JIT compilation. The table also includes 
the execution times using optimizations (opt.) which speeds 
up execution times at the expense of longer compilation 
times. The SSA value denotes the average time for a single 
run while the values in the brackets denote the extrapolated 
value for an ensemble of 50,000 independent realizations 
needed to generate Fig. |4jb). Benchmarks were performed 
on MacOS 10.7, Core 2 Duo @1.4Ghz (64Bit) using one 
core. 



III. Methods 

We consider a general reaction network confined in a 
volume Q under well-mixed conditions and involving the 
interaction of N distinct chemical species via R chemical 
reactions of the type 



sijXi 



■snjXn^ rijXi^ ...^tnjXn, (6) 



where j is the reaction index running from 1 to R, Xi 
denotes chemical species i, kj is the reaction rate of the j^^ 
reaction and Sij and r^j are the stoichiometric coefficients. 
Note that our general formulation does not require all 
reactions to be necessarily elementary. 

The CME gives the time-evolution equation for the prob- 
ability P(n, t) that the system is in a particular mesoscopic 
state ft = (ni, un)^ where is the number of molecules 
of the i^^ species. It is given by: 



dP{n,t) 
di 



N 



1 ]aj {n,n)P{n,t), (7) 



where Sij 



Sij, dj{n^Q) is the propensity function 



such that the probability for the y'^ reaction to occur 
in the time interval [t^t + dt) is given by dj{n^Q)dt 
and E~^'^ is the step operator defined by its ac- 
tion on a general function of molecular populations as 



)i 



The CME is typically intractable for computational pur- 
poses because of the inherent large state space. iNA cir- 
cumvents this problem by approximating the moments of 
probability density solution of the CME systematically using 
van Kampen's SSE |2|, |14j|. The starting point of the 
analysis is van Kampen's ansatz 



n 



[x]^n-^^^e, 



(8) 



by which one separates the instantaneous concentration into 
a deterministic part given by the solution [X] of the macro- 
scopic REs for the reaction scheme ^ and the fluctuations 
around it parametrized by e. fluctuations around it. Note 
that / = lirriQ^oo^/^- The change of variable causes 
the probability distribution of molecular populations P(n, t) 
to be transformed into a new probability distribution of 
fluctuations n(e, t). The time derivative, the step operator 
and the propensity functions are also transformed (see for |j4| 
for details). In particular the propensities expand in powers 
of and are given by 



ttj (n, ri) 



{[X]) 



Note that here we have used the convention that twice 
repeated Greek indices are summed over 1 to A^, which we 
use in what follows as well. Consequently the CME up to 
order rt~^ can be written as 
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where the operators are 
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and the SSE coefficients are given by 
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J-^ ij..r 

d d 



SikSjk---Srk fk^\[^])^ 
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(a) SSE (b) SSA 

Figure 4: Screenshots of iNA's SSE using the EMRE and lOS analysis versus the average of 50, 000 stochastic simulations 
of the gene oscillatory network which exhibits amplified synchronous oscillations which are not captured by the REs (see 
Fig. 8 in Ref. [HI). The EMRE mean concentrations together with the lOS variances in panel (a) are in good agreement with 
SSA simulations in panel (b). The computation of (a) takes less than two seconds compared to more than one hour for the 
stochastic simulation in (b), see Tab. Ill| which highlights the efficiency of the system size method implemented in iNA. 



Note that the above expressions generalize the expansion 
carried out in Ref. |15| to include also nonelementary 
reactions as for instance trimolecular reactions or reactions 
with propensities of the Michaelis-Menten type |16|. Note 
also that the Q^^'^ term vanishes since the macroscopic REs 
are given by dt[Xa] = Efc=i Sakfi^\[^]) leaving us with 
a series expansion of the CME in powers of the inverse 
square root of the volume. In what follows we shall omit 
the upper index in the bracket of the SSE coefficients in the 
case of n = 0. 

Linear Noise Approximation: We now proceed by con- 
structing equations for the moments of the e variables. We 
follow the derivation presented in Ref. f\?\ and expand the 
probability distribution of fluctuations n(e, t) in terms of the 
inverse square root of the volume 



n(6,t) = ^n,(e,t)(^-^-/2, 

J=0 



(15) 



As a consequence there exists an equivalent expansion of 
the moments 

oo 

(e/ee^..e^) = ^[efee/...e^]_^-^]"-^/^, (16) 

i=o 

where the following definition has been used 

[ekei...em]j = J ekei...em^j{e,t)de. (17) 

Using the expansion of the probability density, Eq. ( p3] ), 
together with Eq. (fTOl) and ^TT\ , we find after equating all 
terms of order 



d_ 
dt 



no(6,t) = /:^«)no(6,t) 



(18) 



which is a Fokker-Planck equation with linear drift and 
diffusion coefficients, also called the Linear Noise Approx- 
imation. If the initial state is deterministic, i.e., fi/ft = [X] 
initially, the solution is a multivariate Gaussian distribution 
|2| centered around [e]o = for all times. The time evolution 
equations of the second moment are obtained by multiplying 
the latter by and CiCj respectively and performing the 
integration over e: 



d 



1 



g^[^^^J]0 = J^Mo + ^D,, + (Z ^ j), (19) 

where {i ^ j) denotes the permutations of indices. It then 
follows that within the LNA all higher even moments can 
be expressed in terms of the second moment by virtue of 
Wick's theorem y/7J: 



[eie2...e^]o 



all possible pairings 
P of {1,2. .,C} 



Note that all odd moments are zero since no(e, t) is cen- 
tered. 

EMRE and lOS approximations: The system size expan- 
sion can be used to calculate higher order corrections to the 
moments. The leading order correction to the first moment 
is given by 



1 



^t[e^]l [eji + ^J^'[eae^]o + (21) 



which is equivalent to the EMRE |T8| implemented in the 
previous version of iNA |4 |. The correction to the second 
moment has been derived in detail in Ref. [IS] and is given 



by 



(22) 



which is inherently coupled to the leading order non- 
Gaussian correction of the third moment 



1 



(23) 



Note that the above equations depend on the fourth moment 
which is readily evaluated using Wick's theorem: 

(24) 

In order to relate the above moments back to the moments 



of the concentration variables we use Eqs. ([8]) and to 
find expressions for the mean concentrations and covariance 
of fluctuations which are given by 



}) = [x,]^n-'[e,]i^o{n-^) 



(25) 



(26) 

The term of the latter equation gives the LNA estimate 
for the covariance while including terms to order ft~'^ gives 
the lOS (Inverse Omega Squared) estimate. 



Note that by inspection of Eqs. ( [21] ) and ([22]) it fol- 
lows these corrections are nonzero if the reaction network 
involves at least one bimolecular reaction. Moreover we 
can also estimate the next order correction to the mean 
corrections, i.e., the correction to the EMRE, as has been 
done in Ref. |T5| 

1 



^r^^[e«e/3e^]i + (27) 



3! 



(l)«r 



The solution of the above equation allows us to obtain the 
mean concentration accurate to order 



[Xi] + + n-^[e^]2 + o{n-^), (28) 



where the first term is the RE solution, the first two terms 
constitute the EMRE estimate and including all three terms 
gives the lOS estimate for the mean concentrations. 
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Figure 5: Michaelis-Menten reaction. We have verified the 
soundness of our numerical implementation by comparing 
with the analytical result using the lOS derived in Ref. | [T9| . 
The graph shows the ratio of lOS and LNA variance of 
substrate fluctuations against the ratio S of the free enzyme 
concentration and the total enzyme concentration in steady 
state conditions. This is also compared to the SSA where 
the ratio of SSA and LNA variance has been used. 

Summarizing the system size expansion is obtained by 
defining the vector fssE = (^"^^lna,emre, ^"^^los) with 
components 

^LNA,EMRE = VCC ([Cie^Jo, [ei]i) , 

^los = vec {[eiejekei]o, [eiCjCk]!, [eiej]2, [ei]3) , (29) 

where the symbol vec(-) denotes the vectorization with 
respect to the independent components of the involved fully 
symmetric tensors. Then the vector xlna,emre defines the 
LNA covariance together with the EMRE and xjos defines 
the corrections to the LNA covariance and third moments of 
the SSE together with the order correction to the EMRE 
mean concentrations. It then follows that xsse satisfies the 
linear ODE 

dxssE 



dt 



(30) 



where the upper block triangular matrix B = B{[X]) 
and the full vector A = A([X]) are determined by Eqs. 
( 19|21|22|23|27 ) and depend parametrically on the solution 
of the REs. 

In Fig. [5] we have verified the correctness of iNA's numer- 
ical implementation of the lOS approximation by comparing 
with the compact lOS analytical expression recently derived 
for the case of a simple enzyme catalyzed reaction (see Eq. 
(74) in Ref. [19]). 

IV. Discussion 

In this article we have introduced and implemented the 
lOS approximation in the software package iNA. This allows 
the mean concentrations and variances to be determined 
accurate to order an approximation which is superior 



to the previously implemented methods of LNA (variances 
accurate to order ft^) and EMRE (mean concentrations 
accurate to order As we shown this increased accuracy 
is desirable to accurately account for the effects of intrinsic 
noise in biochemical reaction networks under low molecule 
number conditions. In particular, we have demonstrated 
the utility of the software by analyzing an example of 
gene expression with a functional enzyme and showcased 
the efficiency of the method using a more complex gene 
oscillatory network. We have also extended iNA by a more 
efficient JIT compilation strategy in combination with im- 
proved numerical algorithms which offers high performance 
and enables computations feasible even on desktop PCs. 
This feature is particularly important when analyzing noise 
in reaction networks of intermediate or large size with 
well- separated timescales and some bimolecular reactions, 
conditions that have been shown to amplify the deviations 
from the conventional rate equation description |20|. 

We conclude by pointing out that while stiffness arising 
from multiple timescales in biological networks is still 
an insufficiently resolved problem of stochastic simulation 
methods ||3|, | [2T| , it does not pose significant problems for 
the LNA, EMRE and lOS approximation methods. This is 
since the implementations of the latter in iNA are able to 
deal with stiffness natively using well established implicit 
methods developed for ordinary differential equations. 
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